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SUMMARY 

Wave rotors represent one of the promising technologies for achieving very high core temperatures 
and pressures in future gas turbine engines. Their operation depends upon unsteady gas dynamics and as 
such, their analysis is quite difficult. This report describes a numerical model which has been developed 
to perform such an analysis. Following a brief introduction, a summary of the wave rotor concept is given. 
The governing equations are then presented, along with a summary of the assumptions used to obtain 
them. Next, the numerical integration technique is described. This is an explicit finite volume technique 
based on the method of Roe. The discussion then focuses on the implementation of appropriate boundary 
conditions. Following this, some results are presented which first compare the numerical approximation 
to the governing differential equations and then compare the overall model to an actual wave rotor exper- 
iment. Finally, some concluding remarks are presented concerning the limitations of the simplifying 
assumptions and areas where the model may be improved. 


1.0 INTRODUCTION 

The performance of modern aircraft gas turbines is limited primarily by the pressure ratios and 
turbine inlet temperatures that can be achieved. The wave rotor, when used as a core to a high pressure 
engine, is receiving renewed attention as a potential means of overcoming these limitations using conven- 
tional materials. The wave rotor is a device which utilizes unsteady waves to transfer energy directly to 
and from the working fluid through which the waves travel. As shown schematically in figure 1, it is a 
relatively simple device composed primarily of a series of constant area tubes that rotate about an axis. 
Through rotation, the ends of the tubes are periodically exposed to various ports in which the flow is 
steady, but which initiate or extinguish unsteady waves within the tubes. The ports are circumferentially 
arranged in such a way that those on one end open or close at precisely the time that a wave initiated at 
the other end has traversed the tube. Because of the unsteady nature of the device, each tube of the wave 
rotor is periodically exposed to both hot and cold flow over roughly equal time periods. The mean tem- 
perature of the rotor material may therefore be expected to remain considerably below the peak cycle 
temperature. This characteristic, and the relatively efficient nature of wave compression and expansion 
are the chief advantages of the wave rotor. 

Unlike the conventional gas turbine, which can be successfully analyzed in the steady reference 
frame, the wave rotor is an inherently unsteady device. The flow in the tubes can be approximated by 
the equations of motion for one-dimensional, inviscid, compressible flow of an ideal gas (e.g., the Euler 
equations). These time dependent, hyperbolic equations have closed form solutions only for the simplest 
of flows. Any other flows, such as those found in many practical wave rotor cycles require approximate 
numerical integration. Accurate integration of the Euler equations however, with their discontinuous 
solutions (e.g., shocks and contact discontinuities), is a formidable problem. To this end, the following 
report describes a numerical model which has been implemented at Lewis Research Center expressly for 
the calculation of wave rotor cycles. The model solves the flow equations in a hypothetical wave rotor 
with up to six ports using Roe’s (ref. 1) approximate Riemann solver in an explicit scheme. The ports 
are capable of accommodating supersonic or subsonic inflow or outflow. The model is versatile in that it 
can compute both on and off design flows. 



This report is intended for readers who are at least cursorily familiar with wave rotors, fluid mechanics, 
and elementary fluid computational techniques. Basic discussions on wave rotors may be found in refer- 
ences 2 to 5. For reference however, a brief description of the wave rotor concept will be given in section 
one. Excellent reviews of computational techniques may be found in references 6 and 7. The model 
description will include a presentation of the governing equations and the assumptions used to obtain them, 
a detailed description of the numerical integration technique, a discussion of appropriate boundary conditions, 
and an explanation of the averaging procedure used to obtain the stagnation conditions in the outflow 
ports. Results will then be presented which compare the numerical technique to some known solutions of 
the governing equations. The complete model will then be compared to a simple wave rotor experiment 
(ref. 7). Following the presentation of results, some concluding remarks will be made regarding the valid- 
ity of the simplifying assumptions and some areas where improvements can be made. 


2.0 THE WAVE ROTOR CONCEPT 

The essence of the wave rotor is a series of tubes of constant cross section arranged circumferen- 
tially about an axis to create an annular ring. This is shown schematically in figure 1. The tubes rotate 
on a central shaft and are surrounded by a housing. At the ends of the housing, and perpendicular to the 
tubes, are stationary plates with slots located at various angular positions. The slots, which are the same 
height as the tubes, compose the ports through which flow either enters or leaves the rotor. The primary 
purpose of the rotor rotation is to cyclically expose each tube to the various ports (i.e., timing), thus 
creating the unsteady waves by which the device operates. The easiest way to understand the wave rotor 
process is to adopt a frame of reference on an individual tube and “unwind” the tube from the circumfer- 
ence of the rotor. Each port is then considered to be a reservoir which can be instantaneously opened to 
the tube or closed from it by means of a removable wall (e.g., the stationary end plate). This is some- 
times referred to as the shock tube analogy. With this viewpoint, consider the hypothetical compression 
process shown in figure 2. The gas in the tube is initially at rest at pressure and temperature P 2 and 
T 2 and isolated from both end reservoirs by walls. The left reservoir contains gas at conditions Pj and 
T x which are higher than those in the tube. The right reservoir contains gas at conditions P 3 and T 3 
which are at a particular value to be disclosed momentarily. The wall on the left end of the tube is then 
instantaneously removed. This generates a shock wave which travels to the right at some velocity U g j loc j c . 
Across the shock there is a pressure and temperature rise; however, in general this temperature will not 
match that of the so called “driver” gas that was in the left reservoir. The result of this mismatch is a 
temperature or density discontinuity across which the pressure and velocity are constant, (b). The con- 
tact discontinuity travels to the right at a velocity U which is less than U 8 j loc | c . It is assumed here 
that there is an instantaneous, isentropic acceleration from the resting stagnation state of the tank to the 
static state in the tube when the wall is removed. Strictly speaking, this is not true; however, as pointed 
out by Foa (ref. 2) it is a reasonable first order simplification. The shock wave travels down the tube 
compressing, (i.e., doing work on) the gas in front of it. When the shock reaches the right end of the tube, 
(e), the right hand wall is removed. If the conditions in this reservoir, P 3 in particular, exactly match 
those behind the shock, then the shock will not reflect. The right end of the tube will thus see a steady 
flow until the contact discontinuity arrives some time later, (f). When the contact does arrive, the right 
hand wall is instantaneously replaced. This generates a second shock which travels upstream (left) even 
though the flow itself Is moving to the right, (g). This is often referred to as a hammer shock. Finally, 
when this shock reaches the left hand end of the tube, (i), the wall here is also replaced. Since the gas 
behind the shock is already at rest, replacing this wall does not create a reflection. The gasjn thetubejs 
thus completely at rest in a uniform, compressed state. There are two salient features of this process to 
note. First, the wave speeds are different from the gas speeds. Second, although the processes in the 
tubes are unsteady, as long as the walls are removed the flow at the tube ends (i.e., to or from the ports) 
is steady. A similar series of diagrams may be drawn for the expansion process and these are shown in 
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figure 3. Here, there is a similar gas at rest at state 2. Pj and T l are now lower than the tube. In 
this case there is no right hand reservoir, only a solid wall. When the left wall is removed, an expansion 
wave is initiated which travels to the right, (b); however, the flow itself is to the left. Unlike compression, 
entropy considerations dictate that a discontinuous expansion shock cannot be supported and the wave 
begins to spread. The wave, or fan, is typically represented by a head and tail which represent the front 
and back of the wave. These two parts travel at different speeds through the gas with a continuous change 
in properties between them. In front of the head and behind the tail however, the conditions remain uniform. 
With the wall present at the right end, the expansion will reflect and travel back to the left, (e), spread- 
ing all the while. When the head reaches the left end of the tube, (f), the wall is replaced. Unlike the 
compression scenario; however, the nonuniform properties caused by the spreading of the fan prevent the 
gas in the tube from ever recovering a uniform state. Consequently, from this point in time there will 
always be waves present in the tube. This is an important point concerning expansions and will be dis- 
cussed more fully in subsequent sections. Once a fan is initiated, it is impossible to avoid reflections 
using only instantaneous wall opening and closing. 

The concepts required to understand the wave rotor are now complete. The rotor can be thought 
of as many shock tubes all going through the same process. The opening and closing of walls is accom- 
plished by the ends of the tubes moving from a wall section of the end plates into one of the port regions 
or vice versa. The rotation of the tubes is thus entirely for timing purposes rather than for any aero- 
dynamic effect (actually, some flow turning can be achieved if the tubes are canted but this is not relevant 
to the wave concept). In each port, like the reservoirs described above, the flow is steady. If a glass wave 
rotor could be constructed and the gas could actually be made to change color as, say, temperature or 
pressure change, an observer looking at an operating wave rotor would see a standing wave pattern through- 
out the machine. For a typical wave rotor topping cycle the hot, high pressure gas which drives the 
compression process comes from a burner, the compressed gas behind the shock described above goes to a 
burner, and the expanded gas goes to a turbine. Precompressed air is introduced to the wave rotor via a 
port from an upstream compressor. 

The final concept to be presented, and one which makes the design of wave rotor cycles difficult is 
the requirement of periodicity. Regardless of which ports the tube is exposed to as it rotates on the axis, 
it must end the revolution in the same state at which it began. That is, the cycle must be periodic. This 
allows the port flows to be steady. Consider a gas initially at rest in the tube. Suppose that this gas 
undergoes the compression process described above. In order to complete the rest of this wave cycle it is 
necessary that the gas now compressed in the tube be expanded and scavenged from the tube. Further- 
more, a new charge of gas must enter the tube and be brought to rest at exactly the same conditions at 
which the compression process was initiated. From the reference point of an observer outside the rotor, 
with circumferential velocities neglected, this requirement ensures that the sum of the steady state mass, 
momentum, and energy fluxes through the ports must be zero. From this brief description, it may be 
seen that the wave rotor is a steady flow device (as seen from the outside) in which work exchange is 
obtained by means of unsteady pressure waves. By contrast, most conventional steady flow work 
exchange devices (e.g., turbines and compressors) rely on aerofoils. 


3.0 MODEL DESCRIPTION AND NUMERICAL TECHNIQUE 

3.1 Simplifying Assumptions and Governing Equations 

When discussing modeling of any process there is always debate over which aspects are fundamen- 
tal and must be considered, and which are secondary to understanding and can be neglected. The wave 
rotor is no exception, and the present choice of assumptions represents the author’s opinion on this 
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matter, at this time. The shock tube analogy was preserved for this analysis, and the wave rotor model 
is essentially an integration over time and space of a single tube as it encounters the various ports. This 
simplification implies that no accounting is made for the effects one tube might have on another. It is 
also assumed that the circumferential velocity of the tube is negligible compared with the axial component. 
Thus, calculations do not consider the effects encountered when the velocity of the gas relative to the 
tube is not directly aligned with the tube. The ends of the tube are considered opened to and closed from 
the ports instantaneously. In reality of course, the opening time is determined by the width of the tube 
and the speed of rotation. The flows in the ports are also assumed to adjust instantaneously to those in 
the tube. The argument for this assumption is based on the fact that waves traveling into the ports (which 
arc much larger than the tubes) after the tube opens lose their strength in proportion to the square or the 
cube of the distance that they travel. Hence, in a few tube diameters, or the time corresponding to the 
wave traveling that distance, the strengths of the waves are reduced to mere acoustic perturbations. On 
top of these assumptions are placed the restrictions of one-dimensional, adiabatic, inviscid, perfect gas 
flow. 


The equations that result are usually referred to as the Euler equations and may be presented in 
conservative form as follows: 


— (p) + — (pu) = 0 C 1 ) 

dt dx 

— (pu) + — (p + p u Z ) = 0 (2) 

dt cbc 


i(E) + ±[u(p + E)] = 0, (3) 

dt dx 


E — p 


e + 


2J 


Here p, p, and e are the pressure, density, and internal energy, respectively, 
be used to write the internal energy as, 


e 


£ 


P 

9 


The perfect gas law may 


( 4 ) 


where 7 is the ratio of specific heats. These equations are typically written as a condensed vector 
equation, namely 

— (w) + — (Fl = 0. (5) 

dt ~ dx “ 


Although it is not obvious from equations (1) to (3), it can easily be shown that the flux vector F is a 
function of the conserved vector w, i.e., F = F(w) . 
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3.2 Numerical Scheme 


In the present work, the numerical method used to integrate the Euler equations is an explicit, 
second order, Lax-Wendroff type TVD scheme based on the method of Roe (refs. 1 and 8 to 10). The 
following section of the report will elaborate on this statement. 

If the space of the wave rotor tube (or any space of uniform cross section) is divided into discrete 
“cells” with centers at iAx and faces at (i±l/2)Ax, as shown in figure 4, then Lax (ref. 11) showed 
that the best scheme for integration would be 


n-j-l n 

w. = w. 
— 1 — 1 




At 

Ax 


( 6 ) 


where At is the discrete time step, n is the time index and f is a numerical estimate of the actual flux 
vector F at the cell face (i.e., id=l/2) which is based on information in nearby cells. Although equation (6) 
has the appearance of a first order approximation, it will be seen that the accuracy can be much higher, 
depending on the particular form of f. Equation (6) is consistent not only with the differential form of 
the conservation equation (5) but with the integral form, which may be obtained by integrating equa- 
tion (5) over all x and bringing the time derivative outside the integral to yield 


xmax 

- f (w dx) = -(F -F . )• ( 7 ) 

J — * — xmax —xmin' v 1 

xrain 

This equation truly defines the term conservative for it states that the only change to the amount of any- 
thing (mass, momentum, energy) in the interior of a given space arise from the flux of the conserved vari- 
able and the pressure force at the ends. Unlike equation (5), equation (7) is valid even when w becomes 
nondifferentiable, as in the presence of discontinuities. The solution of equation (7) is often referred to as 
the weak solution. In the discretized system, the same result is found if a summation of equation (6) is 
performed over an interval from i . to L av . Each of the flux terms effectively “telescope” until only 
those at each end of the summation remain. In this line of thinking, it is best to envision the point i as 
representing, not the state of the fluid at that point, but the integral of the state over the width of the 
discrete “cell.” Note that since Lax’s equation (6) satisfies both the differential and integral forms of the 
discrete conservation equations, it too remains valid even in the presence of discontinuities such as shocks. 

The question remains of how to obtain the flux estimates in equation (6) at the cell faces, since the 
only information available lies in the interior of the cells. Two obvious requirements for this flux are that 
the estimate should be based on data in nearby cells, and that if no difference exist in the conserved vari- 
ables of nearby cells (which by definition makes their fluxes identical), then the flux face estimate should 
reduce to the common cell flux. These requirements may be stated as 

$+ 1/2 = ^ 
f(52j. k v..Wj +k ) = E(w) if w. k - Wj +1 . k = ...w. +k 

A less obvious requirement, but one which is of critical importance to the accuracy of the solution, is that 
the estimate should recognize the highly directional nature of the governing equations. Hyperbolic partial 
differential equations propagate information along specific trajectories in space and time, usually referred 
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to as characteristics (refs. 2, 6, 7, 17, and 18). This may be seen mathematically by recasting equation (5) 
in the quasi-linear form, 


Here, the substitution 


dw r , dw 

-= + lti-= = 0. 

dt dx 


dF dF dw dw 

dx dw dx dx 


( 9 ) 


( 10 ) 


has been used. The Jacobian, [A] plays a critical role in the discussion to follow. By definition, the 
Jacobian of any hyperbolic partial differential equation has real eigenvalues and may be diagonalized as 

[A] = [R][A][L], (U) 

In this equation [R] is a matrix whose columns are the right eigenvectors of [A], [A] is the diagonal matrix 
of eigenvalues, and [L] is a matrix whose rows are the left eigenvectors. Furthermore, [Rj is the inverse 
of [L]. Substituting equation (11) into equation (9), and premultiplying by [L] yields 


i dw n i, dw 

l k ^Z + A k l k _z: = 0, 

dt dx 


( 12 ) 


where the l k are the rows of [L] corresponding to each eigenvalue, A k . There are as many eigenvalues 
as there are elements in the conserved vector. If a characteristic trajectory dx/dt — A is followed then 
it is seen that equation (12) becomes a total derivative and the system of partial differential equations is 
reduced to a set of ordinary differential equations along this path, i.e., 


. dw 

l k ^ = 0. 
dt 


(13) 


The eigenvalues are thus speeds at which signals travel in x-t space (for the Euler equations they are 
u+a, u, and u-a where a is the sound speed). By definition of the scalar vector product, the left eigen- 
vectors are directions in conserved space along which the projection of the conservative vector is constant. 
For linear equations, the eigenvectors, and eigenvalues are independent of w. This means that along the 
characteristic paths, which are straight lines in the linear case, a particular linear combination of the con- 
served variables remains constant. If an initial profile of the conserved vector is given, then the profile at 
any other time may be found by simply transporting these linear combinations along the characteristics. 
For any point in the x-t space therefore, a piece of the complete information necessary to define the state 
comes from each of the characteristic paths which reach that point from an earlier time. This is illustrated 
in figure 5, which shows a linear three wave system. For nonlinear equations, the characteristic paths 
and eigenvectors are themselves functions of the conserved variables (and are thus curved), and it is a 
combination of ordinary differential equations which are carried along the path (ref. 2). Nevertheless, it 
is evident that information from one part of the x-t space moves in a specific fashion and has a limited 
range of influence on another part. Any numerical scheme which accurately models the governing equa- 
tions will therefore account for this aspect, which is usually termed “upwinding.” 

Many upwinding methods have been described in the literature; however, the method proposed by 
Godunov (ref. 12) is perhaps the most physically intuitive. For this reason, it is well suited to clarify the 
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upwinding concept. More importantly however, it was the precursor to an entire class of upwinding tech- 
niques, often called Riemann Solvers, of which Roe’s method (which is used in the present wave rotor 
model) is a member. Godunov proposed that each discrete cell be thought of as having a constant value 
of the conserved variable equal to the integral average from i — 1/2 to i + 1/2. As a result of this discre- 
tization, there is a discontinuous jump in conserved variables at each cell face, as shown for a single vari- 
able in figure 6. An initial value problem consisting of two different, semi-infinite, constant states separated 
by a discontinuity is the well known Riemann problem, The solution, which is of the weak form, consists 
of m characteristic signals propagating from the original discontinuity at different speeds. Here m is 
the number of components in the conserved vector (three for the one-dimensional Euler equations). In 
general, the signals consist of a compression shock, an expansion fan, and a contact discontinuity and each 
separates two constant states. This is illustrated in figure 7 which shows the solution of a Riemann prob- 
lem with a shock traveling leftward, a contact, also traveling leftward, and an expansion running right- 
ward. The subscripts 1 and 2 in this figure denote the new constant state regions. For the expansion 
fan, entropy considerations dictate that the resulting signal is not a single characteristic but an infinite 
number of them, representing the continuous change of values from the left to the right state (fig. 3). 
Surprisingly, the solution to the Riemann problem may be found by solving several algebraic relations 
simultaneously (ref. 6). Unfortunately, they are implicit and require an iterative procedure. Nonetheless, 
given the solution to the Riemann problem for every set of adjacent cells, the integral value of the con- 
served variables in a cell at the next time level may be obtained by one of two methods. The first is to 
obtain cell face flux estimates by locating the constant state region (fig. 7) in which the face lies, calcu- 
lating the flux values from the conserved variables there, then using these values in equation (6). This 
method is entirely adequate if the solution to the Riemann problem is exact and if no more than first 
order accuracy is needed. Obtaining the exact solution is computationally expensive however, and it will 
be shown that approximate solutions yield nearly identical results. Furthermore, even with the exact 
solution, first order accuracy is rarely sufficient for practical applications and, as will be seen below, the 
method described above does not retain the information needed to yield better results. Therefore, with 
the solution to the Riemann problem known, consider first the region from iAx to (i+l/2)Ax in figure 7 
at the new time (n+l)At. A region of the space Ax/2 is still occupied by the initial left state Wj. The 
length of this region is Ax/2 + AjAt, where X l denotes the speed of the shock. Similarly, the length of 
the region of space occupied by the state Wj is A 2 At — AjAt, and that occupied by the state w 2 is 
— A 2 At. Thus, the integral of conserved variables at the next time level may be written as 

+ Wj( A 2 At “ *! At ) ~ w 2 (A 2 At). (14) 


(i+irv 


J 


w n+1 dx = w. 


— + A "At 
2 1 


Here, the - superscripts have been added to indicate that these waves have negative speeds. After some 
algebraic rearrangement, equation (14) may be rewritten as 


1 

Ax 


0 + 1 / 2 ) w . 

( w n+1 dx = ! 

J — o 


E 

k=l 


(^k“l^kl) A At 

Ht-’ 

2 k Ax 


(15) 


where, the signal strength vector Aw^ represents the difference between the constant right and left con- 
served states across a given wave (signal). The negative signals have been isolated in using the relation 
A’ — 1/2(A — | A]). Similarly, for a region from i— 1/2 to i, 
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( 16 ) 


1 

Ax 


* w. 

f W n+1 dx = Z! 
J ~ o 

(i-l/ 2 ) 


3 


E 

k=l 


(A k + |A k |) At 

Aw, 

2 k Ax 


Note that the waves and wave strengths in equation (16) pertain to the Riemann problem associated with 
cells i and i-1 while those in equation (15) pertain to i+1 and i. Adding equations (15) and (16) yields 
the desired integral average of w over the cell i at the time level n-H. If one of the signals is an expan- 
sion fan (as in wave number 3 in Fig. 7), then an extra region may be added which has conserved values 
equal to the integral average between the head and the tail of the signal, i.e., 


A^At 

^ 3 ( A 5 - A a) At = J ^ dx - ( 17 ) 

Here, the superscripts h and t refer to the head and tail of the fan, respectively. This average value is 
then incorporated into the solution by adding an extra term to the summation on the right of equation (16). 
With the integral averaged solution obtained in this manner for each cell in the space, the Riemann prob- 
lem may again be solved for every set of adjacent cells, thus allowing the scheme to march forward in time. 

The task remains of insuring that the scheme just described is conservative. To do this, consider 
the Riemann problem illustrated in figure 8. This is the same as figure 7 except two control volumes 
have been sketched in (dashed lines). Equation (7) may be applied to the right hand volume to yield 


— i +1/2 = *, + i 


E A ^ k +A k’ 


(18) 


where the superscript + refers to the positive signals. Similarly for the left hand control volume 


E A ^k A k’ (19) 

Conservation requires that the cell face fluxes of 
equations (18) and (19) be identical. That is, what flows out of the left cell must be exactly what flows 
into the right cell. If this is not the case, the telescoping property, and more formally equation (7), can- 
not hold over the entire space. Equating (18) and (19) results in the requirement 

E i+1 - £i = E A ^ k A k- ( 2 °) 

Furthermore, since equations (18) and (19) are identical, either one, or more easily their average, may be 
used to compute a numerical flux at the cell face which can be used in equation (6), namely 


-i+1/2 -i + 


where the superscript - refers to the negative signals. 


W = 121 

3.2.1 Approximate Riemann solvers and Roe’s method . — The above discussion was prefaced by the 
supposition that the exact solution to the Riemann problem was known for any two adjacent cells. As 
mentioned however, the exact solution is rather computationally costly to obtain and, in light of the fact 
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that the results of the solution are then averaged over the cells, it is not necessary to have such precision. 
The question then arises as to whether an approximate solution to the Riemann problem can be used 
while maintaining the same degree of accuracy in the overall solution. Roe’s method is a particularly 
insightful (though not the only) way of doing just this. 

Roe observed that if the equations of motion are linearized in a particular manner, a weak form of 
the solution can still be maintained, even if the states in two adjacent cells are widely different. In other 
words, his method yields a distribution of, and values for, the constant states in the Riemann problem 
that are different from the true nonlinear solution but have the same integral value. 

Remarkably, the linearization which Roe chooses is derived by assuming small perturbations about 
a mean value of the primitive variables p, />, and u in equation (5). In this linearization, the conserved 
vector and Jacobian have the form 


w = 


* - , ~ * 
p u + pu 

P * +> uu* + ' , ' iJ 


7 - 1 


2 ) 


[A] = 


(7 ~ 3) q2 


1 0 
(3 - 7 )u (7 - 1) 


tl Hu 3 - Hu -(7 - l)u 2 + H 7 u 


( 22 ) 


(23) 


Here, H is the total enthalpy 


*2 -2 

h = _JL_ + Z_, 


7-1 2 

and a is the speed of sound. Note that for linear systems, equations (9) and (5) are identical since the 
Jacobian is constant. It should be kept clearly in mind that, although the linearized equations above are 
derived on the assumption that those quantities with * superscripts are small perturbations of those with 
tildes, that assumption is not made here. The two quantities are to be assumed independent and it is only 
the form of equations (22) and (23) that is used. The linearized Riemann problem can be solved directly 
by integrating equation (9) in time and space in the manner of equation (7). The result is 


i+1 


J w*( n+1 ) dx “ 


*n , *n o 

W. + W. J r 

— i+l — 1 A \~^\ T ~ 7 t *n *Jix 

Ax - £ A k r k l k (w. +1 - w. ) 

1 k=l L 


At, 


(24) 


where the l k and r k are the left and right eigenvectors, respectively, corresponding to the eigenvalues of 
the linearized Jacobian. The summation in equation (24) replaces the term [R][A][L ^w.^ - Wj>. 
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Equation (24) can easily be put in to the form of equation (15) by dividing the integral into two pieces 

centered at i+1/2 and noting that - w*”) = Aw* . Since the system is linear, the possi- 

bility of a fan of any kind does not exist. All the signals take the form of shocks or discontinuities. The 
next step is to replace the linearized conserved vectors on the right hand side of equation (24) with the 

actual conserved vectors from a given Riemann problem. That is, let w. = w. 

By definition, this requires that 


and 


*n 

*i+i 


n 

— M l ' 



Pi 


* 

p i 


= Pi 


♦ Pfii ~ PP 
u. = 

i ~ 


(25) 


JpT n i + yPi + 1 u i + l 

u = 1 — — - 

+ jPi+1 

p — anything, 

where all quantities in equation (25) are evaluated at the time level n. The last three relations in equa- 
tion (25) are found by solving 


/> + u + pu* = pu 


~~ * p* u 2 

/>uu + _ 


pu 

2 


2 


simultaneously at i and i-f-1. The total enthalpy in [A] remains to be obtained before the eigenvectors 
and eigenvalues of equation (24) can be computed. Furthermore, even with this, there is no guarantee that 
the linearized integral on the left hand side of equation (24) approximates the exact solution in any way. 
To achieve this, Roe equates this integral with the exact integral for the Riemann problem over the same 
interval due to the fluxes at i and i+1- That is, 


i+1 

| w( n+1 > dx = 


n . n 
2i+i + 


Ax 


(ll 


i+i -i 


F n At. 


Equating (26) and (24) requires 


*n *n 

• . i — w. 
■l+l —i 


(26) 


(27) 
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which is exactly the relation equation (20) required for the scheme to be conservative. Note that in equa- 
tion (27) the use of the linearized conserved variables is allowed since they are equal to the actual con- 
served variables. The vector relation (eq. (27)) is a set of three algebraic equations which must be solved 

simultaneously; however, there are only two unknowns p and H. Furthermore, it is found that, as in 
equation (27), p remains undefined. Thus, it is not surprising to find that the first two equations of (27) 
are identically solved by the definitions in equation (25), and only the third equation contains new informa- 
tion which yields the needed value of H . Although not necessary in theory, a tremendous simplification 
in the algebra can be obtained if it is required that Au = Au, where the A represents the difference 

between i+1 and i. From equation (25), this is found to be possible if p = • Thus, the change 

in the primitive variables p, p, and u are the same as the change in the linearized primitive variables. 
With this stipulation, the mean enthalpy takes the same form as the mean velocity in equation (25) with 

Uj replaced by Hj. Furthermore, the three scalar values of the products ^(w.^ - w.* n ) , which are a 
measure of the strength of the waves, take on the particularly simple form, 


Ap — paAu 


a 2 A p — Ap 


(28) 




_ Ap + paAu 


2a 2 


corresponding to the three eigenvalues of [A] which are u — a, u, and u + a, respectively. The numerical 
flux estimate to be used in equation (6) becomes 


F. . + F. 


4 + 1/1 - ~2 ' “ \ k ?+ |Akl “ k ' 


(29) 


It is worth noting in passing that since the first two equations in equation (27) yielded the same information 
as equation (25), the resulting equation (eq. (29)) could still be obtained without any assumption as to 
the form of the linearized conserved vector. Only the form of the Jacobian (eq. (23)) is required; how- 
ever, such an approach would not allow the great simplification presented in equation (28). This is one of 
the great strengths of Roe’s approach. Another beneficial feature of Roe’s technique is its ability to cap- 
ture discontinuities accurately. For the particular Riemann problem consisting of a single shock wave or 
contact discontinuity, separating two constant states, and moving with speed S, the Rankine Hugoniot 
relations specify 


Ej. - E, = S(w, - Wj) , (30) 

where subscripts r and 1 denote right and left states, respectively. Comparing this with equation (27) 
requires that S be an eigenvalue of [A] . Thus, if i+1 and i represent two such states, the method 
will single out that signal while setting all the other signal strengths to zero. This amounts to solving 
this particular Riemann problem exactly. 
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3.2.2 Higher order accuracy . — In order to express the ideas in the simplest manner possible, the dis- 
cussion to follow will be restricted to a much simpler differential equation than Euler’s. The advection 
equation, written 


dw f dw 
— + a — = 0, 

dt dx 


(31) 


where a is a constant, is often used as a paradigm for the Euler equations. Although both linear and 
scalar, it is hyperbolic and thus exhibits very similar behavior. Furthermore, the exact solution is simply 
a translation of the initial profile with speed a. This makes comparison with numerical approximations 
particularly easy. Because it is linear, equation (31) cannot develop discontinuities; however, they can be 
prescribed as initial conditions. 

For the advection equation, Roe’s method reduces to 




At 

Ax 


(32) 


where the flux terms have been written out explicitly. In this case, the upwinding concept amounts to 
biasing the space differencing according to the sign of a. To see this, equation (32) may be rewritten as 


n+1 n n n 

w. = w. — a w. — w. , 
1 i\i 1-1 


At ^ n 

a > 0 

Ax 


n-f 1 n 

W. = w. — a 

l i 


In n\ At 

w. M - W. 

' 1+1 1 ' Ax 


a < 0, 


(33) 


In smooth portions of the flow, Roe’s method approximates the true solution of the differential equation 
to within an error margin that is of the same order as the cell size, or time step (Ax, At) multiplied by a 
second order space derivative of the conserved variables. In other words, it is formally First order accurate. 
The accuracy estimate is based on a comparison of the numerical method for cell i with a Taylor series 
expansion about the same point (ref. 6). In the vicinity of discontinuities, all derivatives in space become 
infinite, and the formal meaning of accuracy vanishes. In this case, the accuracy of the solution is in large 
part determined qualitatively by numerical experiment. It is in these cases that Roe’s method has been 
found particularly effective. Most of the x-t space of any problem is occupied by smooth data however, 
and since higher order approximations are possible in smooth regions, the numerical method should recog- 
nize this and attain higher formal accuracy whenever possible. 

All formally second order methods can be derived on the basis of Taylor series expansions about 
the center of the i th cell (ref. 6). Of particular relevance to the present discussion is that of Lax and 
Wendroff (ref. 13). The Lax-Wendroff method may be written as 


n+l 

w. 

1 


£ 

= Wj - - (w i+1 - Wj —1 ) + !L (w i+1 - 2Wj + Wj^j) , 
2 2 


(34) 


where the CFL number v = aAt/Ax has been used. In this and all subsequent equations the super- 
script n will not be used and it is assumed that, unless explicitly stated, all terms refer to the n time 
level. For both equation (32) and equation (34) it is required that \v\ < 1 (in practice, with nonlinear 
problems, this is a strict inequality). This requirement has a formal basis; however, it may be interpreted 
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as insuring that for a given the cell size, the time step is small enough such that a characteristic signal 
does not pass completely through it and deliver information to the wrong cell. The difficulty with the 
Lax-Wendroff method (and all standard second order methods) is that it creates nonphysical oscillations 
in the numerical solution when steep gradients such as shocks are present. To see this, consider that both 
equation (32) and equation (34) can be written in the following form: 

w ” +1 = Wj - C(wj - Wj.j) + D(w i+1 - Wj). (35) 


For the Roe scheme, C = ( 1 / + \u\)/2 and D = (|i/| - i/)/2. For the Lax-Wendroff scheme, 

£' = (*' + * /Z )/2 and D = (i/ 2 — v)j2. For equation (35) to avoid oscillations it is required that 


C > 0, D > 0 and 0 < C + D < 1. 


(36) 


Although these requirements can be proven in a general sense (ref. 14), they are easily illustrated with a 
particular example. Consider the initial distribution of w shown by the solid line in figure 9, with a > 0. 
The exact solution to this problem, as mentioned above, consists of a simple translation of the initial wave 
form with speed a. Since (w-^ — w.) < 0, equation (35) shows that if D is less than zero, the predicted 


value of 


n+1 

w. 

i 


is greater than Wj. Similarly, a simple increment of the indices in equation (35) shows 


that if C is less than 0, then w._^ is predicted less than w i+1 . The numerical result is shown as the 

dashed profile in figure 9 which is obviously oscillatory. From the definitions of C and D, the third 
inequality of equation (36) is simply the limit \u\ < 1. A quick examination of the C and D defined 
above shows that the requirements of equation (36) are met by Roe’s first order scheme, but not by the 
Lax-Wendroff scheme. 


To overcome this problem consider that equations (32) and (34) can be combined by writing equa- 
tion (34) as 


n+l V + \v\ 

W. — W; — 1 1 

1 1 n 


v - V 


( w i “ w i-l) + — [^i+l/2( w i+l - w i) - ^i-l/2( w i - w i-l)] 

(w i+ i-Wj) - [* i+ 1 / 2 (w i+1 - Wi) - < A i _ 1 / 2 (w; - WJ.,)] 


(37) 


The meaning of the <^’s will become clear shortly. Note however, that if <j> = 0, Roe’s scheme is recovered 
exactly and if <j> = 1, the Lax-Wendroff scheme is produced. The subscript on <f> indicates the cell dif- 
ference with which it is associated. The key to making equation (37) successful is to somehow vary <j> in 
such a way that the scheme achieves the highest accuracy possible while still avoiding oscillations. If 
Wj — w._j is factored from the second term on the right of equation (37) and w i+1 — Wj from the third 
term, then the incremental form (eq. (35)) is obtained, with 
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(38) 


C = 


V + \v\ 


1 + 


(1 - y) 


^i+1/2 


r i+l/2 


~ ^i— 1/2 


D 



1 + 


(i - M) 
2 


^i— 1/2 
r i— 1/2 


^i+1/2 f ' 


The quantity r appearing in equation (38) is the dimensionless ratio 


r i+l/2 - 


w i+l-sgn(i/) w i-sgn(i/) 


w i+l 


W: 


(39) 


If a > 0 then D = 0 and the requirements of equation (36) become 0 < C < 1. If a < 0 then C = 0, 
requiring 0 < D < 1. The quantity <f> is dimensionless and must therefore be a function of a dimension- 
less variable. Since the only information available will come from upstream, r is the logical choice, i.e., 


^i+1/2 — ^i-|-l/2( r i+l/2) • 

Application of the conditions in equation (36) to the coefficients in equation (38) requires that 


-2 < ^i+l/2sgn(i/) _ ^ 

1 — |l/| r i+l/2sgn(j/) 


l/2sgn(v) - -T7’ 


which, as may be easily checked, can only be met for any value 0 < \v\ < 1 if 


^i+l/2sgn(v) 

r i+l/2sgn(i/) 


^i-l/2sgn(j/) 


< 2 . 


(41) 


Furthermore, since <f> limits the degree to which the second order scheme is used, it should have a positive 
value. Finally, if r turns out to be negative, which by definition (eq. (39)) means that the solution is 
oscillating, then <j> should be zero in order to recover the Roe scheme completely. With these additional 
constraints laid out, condition (41) can only be absolutely assured if 


M, *(r) < 2, (42) 

r 

where the subscripts have been dropped. These limits ar n shown as the shaded region in figure 10 which 
is a plot of <f> versus r. Also shown is the horizontal lino <j> = 1 representing the Lax-Wendroff scheme. 
Another line, <f> — r is also shown in figure 10. Substitution of <f> = r into equation (38) results in another 
second order (fully upwinded) scheme attributed to Beam and Warming (ref. 15). Recalling that the goal 
of the scheme is to remain second order whenever possible, yet remain oscillation free, reduces the region 
of legitimate <f> to somewhere between the two second order lines but within the restrictions of equa- 
tion (42). This is shown by the hatched lines of figure 10. In the present work, the limiter <f> that is 
used represents the outer limit of the region, namely 
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( 43 ) 


4> = o 

4> = 2r 

<t> = 1 
4> = r 

<j> = 2 

Roe has nicknamed this scheme the “superbee. 
discontinuities. 


if r < 0 

if r < 1/2 

if 1/2 < r < 1 

if 1 < r < 2 

if r > 2. 

' It seems to produce the sharpest resolution of 


The final step is a matter of algebra. Noting that F = aw and that |i/|(l — jj/|) = i/(sgn(i/) — i/), 
equation (37) can be written in the manner of equation (6), with the flux estimate written as 


F:,. + F: 1 1 

f i+l/2 = “ -l a l( w i+i “ w i) + - [sgnM - H^i+i/2 a ( w i+l “ w i)- ( 44 ) 

2 2 2 

The extension of equation (44) to nonlinear systems is fairly straightforward. Each wave system occur- 
ring at the cell face (recall that there are three for Euler’s equation) is considered independent from the 
others. Thus equation (44) assumes the form 


-i+l/2 


1+1 0 1 E I k l A k! a k - \ £ ^lJsgnK) - "kl^V*’ 

2 2 k -_j 2 k=1 


where for each signal, <£ k is evaluated using 


(45) 


a k,i-U/2-sgn(i/ k ) 

r k (46) 

ft k, i+l/2 

The subscripts in equation (46) refer respectively to the wave number and the center of the cell differences 
associated with the A’s in equations (28). Comparison of equation (46) with equation (39) may present 
some confusion since one involves ratios of wave strengths and the other cell differences. It is perfectly 
reasonable however, if it is born in mind that for the scalar advection equation, the solution is simple convec- 
tion of the original wave form. Thus in this case, cell differences and wave strengths are virtually the same. 


3.3 Boundary Conditions 

Up to this point in the discussion, the one-dimensional space in which the equations are solved has 
been considered infinite. That is to say that each point has a surrounding neighborhood from which to 
receive information. In reality of course, there is a finite length to the one-dimensional space and some- 
thing different must be done at the boundaries (for a very complete discussion on this topic see Yee, 
ref. 16). This is clearly illustrated in figure 11 which shows a series of characteristic lines for the advec- 
tion equation (31), with a > 0, in the half space t > 0, 0 < x < x Q . The characteristics shown as 
dashed lines originate from the left boundary and thus require the specification of the solution there, i.e., 
w(0,t) — g(t). At x = x 0 all of the information necessary for the solution is transported there from 
upstream and hence, no boundary specifications are needed. Similar observations may be made for sys- 
tems of linear equations, except only a portion of the information may be specified at the boundaries. 

This is shown in figure 12. Here it is seen that at x = 0, l k w(0,t) = g(t) is specified, where the 
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subscript k corresponds to the positive eigenvalue. The remainder of the information comes from inte- 
rior characteristics running leftward. For nonlinear systems, the situation is similar, but complicated by 
the fact that an ordinary differential equation is carried along characteristics rather than the solution 
itself. Nevertheless, boundaries require that some, but not all of the flow information be specified, and 
what that information is depends upon the direction of the characteristics. For the Euler equations there 
are several possibilities. The discussion to follow will concentrate only on a left hand boundary with the 
understanding that directly analogous arguments may be made for the right hand side. 

Consider the half space where at t = 0, the conserved vector w is at some uniform state for 
x > 0. For t > 0, what can and/or must be specified about the boundary? Either the conserved vari- 
ables or the so called primitive variables p, p , and u (of which the conserved variables are composed) 
may be considered. This discussion will consider the latter. In either case, the answer is that what is 
specified depends upon the outcome, and therein lies a major difficulty to the application of boundary 
conditions. Specifically, the boundary specifications depend on whether the resulting flow is in the posi- 
tive or negative direction, and subsonic or supersonic. To see this, consider the illustration of figure 13 
which shows the possible outcomes when boundary conditions are applied. If the resulting flow is inward, 
there are three possibilities. The first, (a), is that all of the characteristics have positive speeds (eigen- 
values). This is the supersonic condition and requires that all three primitive variables must be specified 
at x = 0 because none of the information from the interior can reach the boundary. The second, (b), 
and third, (c), situations are subsonic inflow conditions which in general create a shock and a contact 
discontinuity, or an expansion and contact discontinuity, respectively. Here, since two of the characteris- 
tics have positive speeds, two primitive variables must be supplied at x — 0. If p and u are specified 
then there are two jump conditions which must be met across the shock or expansion since these variables 
are constant across the contact discontinuity. This is impossible to satisfy and thus it is either p and 
/>, or u and p which must be specified at x = 0. Note that the complete solution for x = 0 in these 
cases relies on knowing the conditions in region 1 (the first cell of the computing domain) as well. This 
constitutes the information that is carried along the left running characteristic from the interior. If the 
flow is outward, there are also three possible scenarios. In the first two, (d) and (e), only one characteris- 
tic has a positive speed and therefore only one boundary variable may be supplied. This is usually the 
pressure. The remaining information is obtained by satisfying the jump conditions across the shock or 
expansion fan that has resulted. The third situation, (f), is that where the tail of the resulting expansion 
fan has zero or negative speed. Here, the velocity at the exit is sonic and nothing can be imposed from 
the outside. The value of the primitive variables at x = 0 are supplied by interior information and the 
condition that u = a. Actually, there is a fourth outflow situation that may occur if the initial flow is 
supersonic. Here again, the outside world has no effect and all information comes from the interior. From 
the above description, the boundary specification difficulty is clearly seen. The boundary conditions, which 
determine the resulting flow, cannot be applied until the resulting flow is known. Some mechanism is 
apparently needed whereby the outcome may be known before the conditions are applied, thus assuring 
that what is specified is correct (i.e., well posed). The author’s present work on boundary conditions for 
the Euler equations, which will be described below, has centered on doing just this. 

For numerical approximations, and those utilizing Riemann solvers in particular, the application of 
boundary conditions amounts to adding a so called “image cell” to the left of the first cell (or to the right 
of the last cell) in the computing domain, as shown in figure 14. At each new time step the correct value 
of the primitive variables at x = 0 is found using appropriate specified data and that from the first cell. 
These values are then used to compute conserved variables which are in turn assigned to the image cell. 
When the subsequent Riemann problem is solved between the image cell and the first cell, the signals 
traveling into the computing space carry the correct information. Those signals directed out of the com- 
puting space (into the image cell) are ignored. 
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The first step in the process is to assume that the flow is outward and assign to the boundary a pres- 
sure. If the pressure ratio of the boundary pressure to that in the first cell, PR = Pq/Pj is less than 1.0 the 
assumption is made that the resulting flow contains an isentropic expansion fan. For one-dimensional 
isentropic flow of a perfect gas it may be shown (ref. 17) that along left running characteristics the quan- 
tity (2/7 — l)a — u is constant. That is 

2 2 (A7\ 

a o “ % = a l “ u l> l 47 ) 

7-1 7-1 

where the subscripts 1 and 0 refer to the first cell and the point x = 0 (or image cell), respectively. 
Noting that for an isentropic, perfect gas 


— = PR 7-1 / 2 ' 1 ', 

a l 


equation (47) may be rearranged as 


Mq = 



M, 


7 - 1 


PR-(7-1/2i) ? 


(48) 


* s 

where M is the Mach number. This equation may be used to obtain two critical values, Mj and 
Mj is the value of Mj which, for the specified pressure ratio, results in Mq — 0. M^ is the value of Mj 
for which Mq = — 1 (i.e., sonic). These values may be written as 


M* = — - (l - PR 7 " 1 / 27 ) , = M* — PR 7-1 / 27 . 

7 - 1 


(49) 


If PR is greater than 1, then it is assumed that the resulting flow contains a shock wave and the normal 
shock relations (ref. 18) may be used to obtain the Mj and another value, M as . M as is the Mach num- 
ber of the flow in front of the shock as measured from the reference frame of the shock itself; however, it 
is also the value of M 1 which yields a standing shock at the exit for the specified pressure ratio (i.e., the 
characteristic velocity is zero). These quantities are 



7 + . 1 (PR - 1) + 1 , M* 

27 


2 

7 + 1 


/ \ 



(50) 


Equations (49) and (50), with 7 ~ 1.4, are plotted in figure 15 for various values of PR. If PR and 
M 1 for a given boundary are plotted on this figure, then the nature of the resulting flow can be predicted 
from the location of the point. If M 1 < Mj for any value of PR then the resulting flow will be outward 
(i.e., negative). If M^ < M x < M 1? the flow is subsonic and only the boundary static pressure may be 
specified. The values of the other primitive variables are found using either equation (48) and the isen- 
tropic perfect gas law, or the normal shock relations, depending on the value of PR. If PR > 1 and 
M x < M aa , or PR < 1 and Mj < - 1 , then the outflow will be supersonic and all primitive variables 
must be obtained from the interior. In the model this is done by simply assigning the value of the first 
cell to the image cell. However it is obtained, it does not affect the computation in the interior of the 
domain. If PR < 1 and — 1 < Mj < M®, then the flow at the exit is exactly sonic and the primitive 
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variables at the boundary are again obtained by equation (48), the isentropic perfect gas law, and the 
knowledge that = — 1. 

If Mj > Mj, then the flow is inward and additional tests must be applied to determine the appro- 
priate boundary specifications. For inflow conditions, either the static pressure and density, or the pres- 
sure density and velocity must be specified depending on the inflow Mach number. This is mathematically 
correct; however, in the laboratory, these variables cannot be controlled. The analogy being used for the 
wave rotor is that of a shock tube connected to reservoirs at each end. These reservoirs are considered to 
contain stationary gas at a specified pressure and temperature (or density) which, by definition, are stag- 
nation conditions. If the flow is from the tube into the reservoir, as in a jet, then the static pressure of 
the jet (for subsonic flow) is the reservoir stagnation pressure and the reservoir temperature is irrelevant. 
If the flow is from the reservoir into the tube however, then the static conditions do not match the stag- 
nation conditions, and it is only the stagnation state that can be controlled. In the present work, this dif- 
ficulty is skirted by assuming that the stationary reservoir gas accelerates instantaneously, and isentropic- 
ally to the static inlet conditions whenever inflow is predicted. A justification for this assumption was 
provided in the introduction and will not be repeated here. For steady, isentropic flow, the relationship 
between stagnation conditions and the Mach number is 


/ /> 
Po 

7 - 1/7 

fa') 

a 0 
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a 0 
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(51) 


where the prime symbol denotes the stagnation state. 


This equation may be rearranged to obtain 
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(52) 


Equations (48) and (50) may also be rearranged to the form 


Uj = _ = Mj + 


7 - 1 


PR 7 


- l] if 


PR < 1.0 


M! - 


7 + 1 
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m! s 


if PR > 1.0, 


(53) 


where Mj SS is defined in equation (50). Equations (52) and (53) may be solved simultaneously for PR; 
however, it is not a closed form solution and some type of iteration must be used. If the solution yields 
an inlet Mach number which is less than 1, then the pro ess is complete. If the calculated Mach number 
is greater than 1 however, then the results may not be used. In this case, the inflow is supersonic and all 
boundary conditions must come from the reservoir. This requires an extra piece of boundary information 
which is the inlet Mach number. With the Mach number and stagnation conditions known, all of the 
required primitive variables at the boundary may be calculated. 


The above discussion has covered boundary conditions for inflow to or outflow from the computing 
space; however, there is a third set of boundary conditions which must be dealt with and these arise in 
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the presence of a solid wall. Here the flow is abruptly brought to a halt. In this case one of the primitive 
variables is already known, namely u Q — 0. The others may be found by noting that in the case of a solid 
wall, M x = M 1 in equations (49) and (50). Which equation is used depends on whether Mj is positive or 
negative, respectively. These equations can then be rearranged to solve for PR and appropriate relations 
(isentropic laws or shock laws) may be used to find /> Q . 


3.4 Port Flows 

The discussion thus far has focused on the gasdynamics of a single tube in the wave rotor as it is 
exposed to the various ports on the periphery. It was pointed out in the introduction however, that all of 
the tubes (fig. 1) undergo the same process. Thus, from the point of view of the ports, the flow is steady. 
Although the internal dynamics are the heart of the wave rotor operation, it is ultimately these port flows 
that are of interest for they are the measurable performance indicator. Specifically, it is the stagnation 
conditions at the various ports that are desired. For ports where the flow is into the rotor, these condi- 
tions are already known due to their being required for proper boundary condition specifications. For 
outflow ports however, all that is specified is the static pressure and the stagnation conditions must be 
estimated by some additional analysis. 

As was the case for the boundary conditions, the discussion to follow will concentrate on the left 
side of the wave rotor. Consider the steady state control volume shown in figure 16. The profile sketched 
on the right side of this figure represents a possible time trace for a flux variable of one of the wave rotor 
tubes as it traverses the port. This is also the steady (but spatially varying) profile observed from the 
port frame of reference. The profile has been deliberately made to look nonuniform in order to highlight 
the generality of the analysis. The profile on the left represents the same flux variable after all the non- 
uniformities have mixed out. This is assumed to occur some distance down the port duct. The particular 
distance is unimportant since, in the present analysis, the effects of wall friction have been neglected. In 
reality, the frictional losses are probably not negligible and a more detailed analysis will be necessary. It 
is reasonable however, to assume that the mixing losses are at least dominant. Assuming that the flow is 
uniform from the top to the bottom of the duct, and that the cross section of the duct remains constant, 
the following steady state mass, momentum, and energy equations may be written: 
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(54) 


Here, T' denotes the stagnation temperature, the overbar refers to the mixed state, the subscript 0 refers 
to the state at the boundary of the tube, and 6 X and # 2 designate the location of the port. Since the 
integrals F masg , F mom , and F en can all be numerically evaluated from the model, they are known and 

equations (53) may be solved simultaneously for p, />, and u using the additional perfect gas relation 
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(55) 
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With the mixed static quantities known the stagnation pressure can easily be calculated. 


4.0 RESULTS 

There are two major aspects of the wave rotor model which must be examined in this section of the 
report. The first is the degree to which the numerical method matches the differential equations which it 
approximates. The second is the validity of the simplifying assumptions on which the equations are based. 
Relatively speaking, the first aspect is easier to examine than the second and it shall be considered first. 

Figure 17 presents six shock tube test scenarios which have known analytical solutions. For review, 
it is noted here that the solutions are analytical under the assumptions of inviscid, one-dimensional, per- 
fect gas flow in a tube with ends that adjust instantaneously to the stagnation conditions present in the 
attached reservoirs. In each case, the exact (analytic) solution is represented by solid lines, while the 
various symbols represent the computed results. The tests each consist of setting the gas in the tube at 
some initial state and then applying boundary condition at t = 0. The results presented represent an 
instantaneous “snapshot” of the state in the tube at some later time. In the figure, the primitive vari- 
ables p, />, and u have been nondimensionalized by arbitrary values p 1? p v and a r Similarly, time 
has been nondimensionalized as taj/L and distance as x/L, where L is the length of the tube. The 
CFL number is defined as Ataj/Ax. In all cases, the gas has a ratio of specific heats, 7 equal to 1.4. 

The tube contains 200 cells; however, only the value of every fourth one is shown in order to avoid 
clutter. The various scenarios will be described below; however, a cursory glance at the six plots clearly 
shows the excellent agreement between analytical and numerical results. As promised, the resolution is 
very high, with no spurious oscillations and the calculated speed of the various waves is nearly exact. 

The first scenario, in figure 17(a), shows the results of a tube initially containing high pressure, 
high density gas in the left half, and low pressure, low density gas in the right half. The two sides are 
initially at rest with a diaphragm separating them. The diaphragm is instantaneously removed a t = 0 . 
This creates a shock and a contact discontinuity traveling rightward and an expansion fan traveling to 
the left. This test does not involve the imposition of any boundary conditions and is therefore a test of 
the numerical method for the interior points alone. 

Figure 17(b) to (d) are the boundary scenarios shown in figure 13. Figure 17(b) is a test for super- 
sonic inflow on the left and subsonic outflow on the right. Here, the gets in the tube is initially in a uniform 
state throughout the tube and traveling to the right subsonically. At t = 0 the right end is exposed to a 
high pressure reservoir which generates a left running shock but does not stop the outflow. The left end 
is also exposed to a high pressure reservoir which generates a right running shock and contact discontinu- 
ity. The reservoir pressure is high enough in this case to create supersonic inflow. Figure 17(c) is also a 
left inflow, right outflow test except that in this case, the gas is Initially at rest, the inflow is subsonic 
and the reservoir on the right is at a low pressure which generates an expansion fan instead of a shock 
running leftward. Figure 17(d) is yet another left inflow, right outflow scenario. Here, the gas in the 
tube is initially moving to the right. At t = 0 the left end is exposed to a reservoir with a lower stagna- 
tion pressure than the gas initially in the tube. The pressure is not low enough to stop or reverse the 
flow and the result is an expansion fan and contact discontinuity combination traveling rightward. On 
the right end, the tube is suddenly exposed to a reservoir at extremely low pressure. This produces an 
expansion fan traveling leftward with a sonic condition at the exit. 
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Figure 17(e) is a test of the boundary conditions imposed by the presence of a wall. The flow is 
initially uniform and traveling rightward. At t — 0 walls are imposed on both ends. This generates a 
left running shock on the right end and a right running expansion fan on the left. 

Finally, figure 17(f) is a test for wave interactions and flow reversals. An x-t diagram accompa- 
nies this figure in order to clarify the process. The gas in the tube is initially at rest with a contact dis- 
continuity positioned at x/L - 0.695. At t = 0 both ends are exposed to high pressure reservoirs which 
generate shocks traveling toward the center of the tube. The stagnation temperatures in the reservoirs 
have been chosen such that no contact discontinuities are created. The two shocks collide precisely at the 
initial contact discontinuity, which results in two new shocks traveling outward and the contact moving 
to the left. When the left running shock reaches the left boundary, the pressure and temperature behind 
it are enough to cause outflow from that end. The result is an expansion fan traveling rightward. The 
plot shows the state of the gas in the tube after the fan has been generated but before any further 
interactions occur. 

It is apparent from figure 17 that, for the assumptions made, the numerical model is an excellent 
approximation to the solution of the resulting differential equations. The question still remains however, 
of whether the assumptions made are valid. That is to say, is the perfect gas, shock tube analogy a valid 
wave rotor model? The easiest way to answer this question would be to compare the model to some experi- 
mental data. Unfortunately, this type of information is extremely difficult to obtain. It seems that those 
experiments which were successful either remain proprietary, were scantily documented, or contained addi- 
tional devices such as burners whose effects on the rotor were inseparable from the rotor performance itself. 
The only useful data which the author has found comes from an experiment conducted by Kentfield (ref. 19). 
He designed and built what he called a “dynamic pressure exchanger.” This wave rotor device split a sin- 
gle incoming flow into two outgoing flows. One stream emerged at a higher stagnation pressure than the 
incoming stream, and one at a lower value. The transfer of energy from one part of the stream to another 
was accomplished solely with pressure waves. The relative mass flow of each of the emerging streams 
depended on the amount of energy transferred. An conceptual view of the cycle is shown on an x-t dia- 
gram in figure 18. Air enters through the port labelled M (medium pressure). This generates a shock 
wave which travels leftward, compressing the stationary gas in the tube and causing leftward flow. When 
the shock reaches the left side, the high pressure port is opened. This creates a right running shock, but 
does not stop the flow from exiting the high pressure port (fig. 17(b)-right side). When this shock 
reaches the right side of the tube, the inflow port is closed. This, of course, creates a left running expan- 
sion (fig. 17(e)-left side) but it is assumed weak and therefore not shown. Next, the low pressure port is 
opened. This creates a left running expansion which is shown as a single dashed line in the figure. When 
this fan reaches the left side, the high pressure port is closed. At this point the wave pattern becomes 
quite complex because the port closing generates a right running shock which is moving through a reflected 
expansion. The shock is also assumed weak and it is the expansion which dominates the field. When this 
expansion reaches the right side, the low pressure port is closed and the cycle is presumably complete 
(even though this too will create a shock), due to the large number of waves which were not accounted 
for. As expected, the experimental performance of this device was poor. Nevertheless, there was a docu- 
mented experiment and figure 18 shows the results. In the experiment, the rotor speed and port locations 
were set and the port pressures were varied over a range of values. The mass flow from the ports was 
then measured. Figure 18 is a plot of the ratio of the high pressure outflow port to the inflow port stag- 
nation pressure versus the ratio of the low pressure outflow port to the inflow port stagnation pressure for 
various families of constant high pressure outflow port to inflow port mass flow ratios. The dashed lines 
represent the experimental results. Also shown in the figure are the results using the wave rotor model. 
The solid circles represent the actual model output. The constant mass flow ratio lines (dotted) were 
interpolated from these results. It is seen that for a given mass flow ratio, an increase in p^ results in a 
lower p£ indicating greater work transfer. For values of Pjj/Pm below 1.2 the model appears to provide 
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reasonable performance estimates. It tends to overpredict performance however, with the error growing 
worse as pressure ratios are increased. This result is entirely expected since the model is an idealization. 
There are several potential reasons for this discrepancy, all of which will be accounted for in future imple- 
mentations of the model. One possibility is the effect of finite port opening and closing times. Each tube 
in the wave rotor has a finite width and does not open or close instantaneously when it rotates into or 
out of the port region. This results in a high loss, jet like flow into the tube until the tube is fully exposed 
to or removed from the port. Another source of error in the model is the lack of an accounting of viscous 
effects. These would tend to degrade the performance in the manner exhibited by the experiment. A third 
error source of disagreement is the effect of leakage flows in the device. A finite gap exists between the 
ends of the tubes and the face plates of the wave rotor. High pressure gas in the tubes leaks through this 
opening to the low pressure region outside the rotor. In all probability, the disparities between the model 
and the experiment are caused by a combination of all three of these error sources. F urthermore, there 
are several additional error sources which have not been mentioned and will not be discussed (e.g., duct 
losses, circumferential velocity effects, etc.). 

The model data for figure 18 was somewhat difficult to obtain and warrants some further discus- 
sion. In order to calculate any cycle, an initial state is needed for the gas in the tube. It was not possible 
to obtain this from the experimental results and therefore, a rather “brute force” technique was adopted. 
The initial state was simply chosen arbitrarily. In particular it was set as uniform and equal to the state 
of the inlet port reservoir. The tube was then run past the port boundary conditions repeatedly until the 
state of the gas in the tube at any point in time exactly matched the state one cycle earlier. This amounts 
to obtaining a fixed wave pattern in x-t space. It was found that 15 cycles of the tube were sufficient to 
insure this condition. One such pattern is shown in figure 19. This is an x-t plot of pressure density 
and velocity using various shades of grey to indicate relative values. This pattern represents one of the 
data points in figure 18. Note the presence of many additional waves besides those shown in the insert of 
figure 18. 


5.0 CONCLUDING REMARKS 

A model describing the dynamic performance of the wave rotor has been developed and was 
described above. This model is based on an analogy between the tubes of the wave rotor and the sta- 
tionary shock tube. The differential equations which dominate the behavior of the gas flow in such tubes 
are the one-dimensional Euler equations. It was found that the numerical techniques used in the model 
provided an excellent approximation to these equations. Furthermore, a comparison with the limited 
experimental wave rotor data available indicated that the shock tube analogy is reasonable. It was found 
however, that some of the simplifications associated with the analogy need to be accounted for. In par- 
ticular, it appears that frictional effects are more significant than was originally believed. These may be 
accounted for with the addition of an appropriate source term to the original Euler equations. However, 
this addition adds a level of complexity to the equations since they are then no longer conservative. Another 
assumption which must be investigated is that of instantaneous opening of tubes. From the literature 
available on past wave rotor experiments it is known that the tubes do not open in this fashion and the 
losses associated with finite opening time appear to be significant. 


REFERENCES - - - ... 

1. Roe, P.L.: “Characteristic Based Schemes for the Euler Equations,” Annual Review of Fluid Mechan- 

ics, 1986, vol. 18, pp. 337-365. 


22 



2. Foa, J.V.: Elements of Flight Propulsion , John Wiley and Sons, 1960. 


3. Taussig, R.T.: “Wave Rotor Turbofan Engines for Aircraft,” Presented at The ASME Winter Annual 

Meeting, New Orleans, LA, 1984. 

4. Croes, N.: “The Principle of the Pressure-Wave Machine as Used for Charging Diesel Engines,” Pro- 

ceedings of the 11th International Symposium on Shock Tubes and Waves , Seattle, WA, 1977. 

5. Weatherston, R.C.; and Hertzberg, A.: “The Energy Exchanger, a New Concept for High-Efficiency 

Gas Turbine Cycles,” Transactions of the ASME , No. 66-GT-117, 1966. 

6. Hirsch, C.: Numerical Computation of Internal and External Flows , John Wiley and Sons, 1990. 

7. Anderson, D.A., et al.: Computational Fluid Mechanics and Heat Transfer , McGraw-Hill, 1984. 

8. Roe, P.L.; and Pike, J.: “Efficient Construction and Utilization of Approximate Riemann Solutions,” 

Computing Methods in Applied Sciences and Engineering , 1984, vol. 6, pp. 499-518. 

9. Roe, P.L.: “The Use of the Riemann Problem in Finite Difference Schemes,” Proceedings of the 

Seventh International Conference on Numerical Methods in Fluid Dynamics , 1980, pp. 354-359. 

10. Roe, P.L., ed.; Numerical Methods in Aeronautical Fluid Dynamics , Academic Press, 1982. 

11. Lax, P.D.: “Weak Solutions of Nonlinear Hyperbolic Equations and Their Numerical Computation,” 

Communications in Pure and Applied Mathematics , 1954, vol. 7, pp. 159-193. 

12. Godunov, S.K.: “A Difference Method for the Numerical Computation of Discontinuous Solutions of 

Hydrodynamic Equations,” Mat Sbornik , 1959, vol. 47, pp. 271-306. 

13. Lax, P.D.; and Wendroff, B.: “Difference Schemes for Hyperbolic Equations With High Order of 

Accuracy,” Communications in Pure and Applied Mathematics , 1964, vol. 17, pp. 381-398. 

14. Sweby, P.K.: “High Resolution Schemes Using Flux Limiters for Hyperbolic Conservation Laws,” 

SIAM Journal of Numerical Analysis , 1984, vol. 5, pp. 995-1011. 

15. Warming, R.F.; and Beam, R.M.: “Upwind Second Order Difference Schemes and Applications in 

Aerodynamics,” AIAA Journal , 1976, vol. 14, pp. 1241-1249. 

16. Yee, H.C.: “Numerical Approximations of Boundary Conditions With Applications to Inviscid Equa- 

tions of Gas Dynamics,” NASA TM-81265, 1981. 

17. Anderson, J.D.: Modern Compressible Flow , McGraw-Hill, 1982. 

18. Thompson, P.A.: Compressible Fluid Dynamics , Department of Mechanical Engineering, Rensselaer 

Polytechnic Institute, 1988. 

19. Kentfield, J.A.C.: “An Examination of the Performance of Pressure Exchanger Equalizers and 

Dividers,” Ph.D. Thesis, University of London, 1963 (also in Journal of Basic Engineering, Sept. 
1969). 


23 



Figure I Wave Rotor Schematic 
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Figure 4 Numerical Cell Nomenclature 
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Figure 5 Characteristic Paths In X-T Space 
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Figure 6 Godunov's Integral Cell Average 
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Figure 1 7 Test Cases 
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